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"It is completely idle to say that the tides are caused by the conflict arising 
between the motion of the earth and the motion of the lunar sphere, not only 
because it is neither obvious nor has it been explained how this must follow, but 
because its glaring falsity is revealed by the rotation of the earth being not con- 
trary to the motion of the moon, but in the same direction. Thus everything that 
has been previously conjectured by others seems to me completely invalid. But 
among all the great men who have philosophized about this remarkable effects, 

I am more astonished at Kepler than at any other. Despite his open and acute 
mind, and though he has at his fingertips the motions attributed to the earth, he 
has nevertheless lent his ear and his assent to the moon's dominion over waters, 
to occult properties, and to such puerilities” 

Galileo Galilei, 

Dialogue concerning the two chief world systems, p. 462. Translation by 
S. Drake. Univ. of California press. 1953. 



CONTRIBUTION TO THE THEORY OF 
TIDAL OSCILLATIONS OF AN ELASTIC EARTH. 
EXTERNAL TIDAL POTENTIAL 


Peter Musen 


ABSTRACT 


In the present article we establish the differential equations 
of the tidal oscillations of the Earth under the assumption 
that the interior of the Earth is laterally inhomogeneous. 

We develop the theory using vectorial and dyadic symbol- 
ism to shorten the exposition and to reduce the differential 
equations to a symmetric form convenient for programming 
and for numerical integration. The formation of tidal 
,r buldges" on the surfaces of discontinuity and the changes 
in the internal density produce small periodic variations 
in the exterior geopotential which are reflected in the 
motion of artificial satellites. The analogues of Love 
elastic parameters in the expansion of exterior tidal 
potential reflect the asymmetric and inhomogeneous struc- 
ture of the interior of the Earth. 
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Introduction 

In the present article we establish the differential equations of the tidal 
oscillations of the Earth under the assumption that the interior of the Earth is 
laterally inhomogeneous . We develop the theory using vectorial and dyadic 
symbolism, partly to shorten the exposition and, partly, to reduce the differ- 
ential equations to a symmetric form convenient for programming and numer- 
ical integration. 

It has been customary in the investigations of the oscillations of the Earth 
to assume the spherical symmetry of internal density and of Lame‘S elastic 
parameters. In addition, the undisturbed Earth has been considered to be in 
hydrostatic equilibrium. However, we have already reached the stage where 
the hypotesis of spherical symmetry cannot explain the results of modern 
seismic observations. The analysis of the seismic information provided by the 
Peru-Bolivian border earthquake (Dziewonski, 1970), shows that the measure- 
ments of the group velocities along the world circling path indicate the existence 
of lateral inhomogeneities in the crust and mantle. Dziewonski also finds that 
the measurement of the surface wave dispersion provides the means for the 

determination of the internal density distribution. From the analysis of the 

!> 

information supplied by the Alaska earthquake Toksoz and Anderson (1966), 
conclude that the available phase velocity data are sufficiently accurate to 
indicate the regional variations and the heterogeneity of the upper 400 km 


mantle. 
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Arkani-Hamed (1970), suggests the lateral variations of density in the upper 
mantle are approx. 3% and of the shear modulus approx. 20%. It is of importance 
to extend our understanding of the mechanism of the Earth's tidal oscillations, 
partly for the sake of the interesting geophysical problem itself and partly 
because the tidal oscillations produce small periodic changes in the geopotential 
which are reflected in the motion of artificial satellites. The problem of tidal 
perturbations in the motion of artificial satellites recently attracted the atten- 
tion of a number of authors (Kozai, 1965), (Newton, 1968), (Kaula, 1969), (Smith 
and al. 1971, 1973), (Musen and Estes, 1972), (Musen and Felsentreger, 1973), 
(Douglas and al., 1972, 1973), (Lambeck and Cazenave, 1973), (Musen, 1973). 

The periodic variations in the geopotential are caused by changes in the internal 
density and by the formation of tidal "buldges" on the surfaces of discontinuity 
inside the Earth and on the free surface (Biot, 1965). Both changes are small, 
but their totality produces appreciable and observable perturbations in the 
motion of artificial satellites. A periodic expansion of the tidal perturbations 
in the height of satellite GEOS-A above the Earth, in terms of the arguments 
of lunar theory, was obtained by R. Estes using the Musen and Estes semi- 
analytical theory of tidal effects (1972) in the motion of artificial satellites. 

They are quite observable. Estes' expansion contains two long period terms 
with amplitudes of 3.8 m. and 5.0 m., respectively. 
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A proper theory of the Earth's tidal oscillations, bodily and surface, 
can help us to understand some characteristics of satellite motion. 

On the other hand, a proper theory of the tidal perturbations of satellites can 
give us some insight into the Earth's internal structure and its elastic properties. 
The average elastic parameters of the Earth for different satellites were re- 
cently obtained by Kozai (1965), Newton (1968), Smith and al. (1971, 1973) and 
Douglas and al. (1972). Typically the same value of the elastic factor is used 
for all harmonics of a given degree in the expansion of the exterior tidal 
potential, irrespectively of the tidal frequency. However, the work of Alterman 
and al. (1959) indicates that the elastic response of the Earth is different for 
different tidal frequencies. More computational work is required in order to 
understand the full implication of this idea on the theory of motion of artificial 
satellites. The parameter of the Earth’s elastic response (Love number) to the 
tidal force of a given frequency, as it appears in the expansion of the exterior 
tidal potential, can be obtained as a by-product of the numerical integration of 
the partial differential equation controlling the tidal oscillations. The integra- 
tion should be extended over the whole Earth. 

In the present article we obtain, starting from the principle of D'Alembert, 
a modified form of the Biot (1965) differential equation for the elastic oscilla- 
tions, using the apparatus of vector and dyadic analysis. We find the form of 
the perturbative terms which disappear together with the initial deviatory stress. 
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The method of Biot includes the effect of the initial stress, as well as the effect 
of the additional stress as created by the vorticity of the displacement field. 

Biot's formulation is Lagrangian. The displacement of a material point is 
given relative to its "initial" position. We incorporate into the differential 
equation of the tidal oscillations of the solid Earth the effect of the semi-diurnal 
nutation and of the geos trophic force (Molodensky, 1953), (Melchior and Georis, 
1968), (Melchior, 1971), and also the effect of the initial deviatory stress. It is 
of interest to note that in the Biot formulation the density factor, associated 
with the forces of inertia, retains its initial value and there is no necessity to 
consider its time variation. The geostrophic force in tidal theory is more 
important than in the theory of free oscillations of the Earth, where it is only 
a perturbative term. The rotation of the Earth causes only a "fine splitting" 
of the spectrum of free oscillations’ frequencies and the analytical perturba- 
tions technique, similar to one in quantum-mechanics, can be applied (Madariaga, 
1972), (Dahlen, 1968, 1972) (Luh, 1973). 

In the tidal theory the forced oscillations take place and a large number of 
significant tidal frequencies are of the same order of magnitude as the angular 
speed of the Earth's rotation itself. As a consequence the geostrophic term in 
the differential equation for tidal oscillations is no longer a perturbative term. 

In our exposition we remain in the frame of ideas of classical mechanics. Like 
Biot, we make use of the principle of D’Alembert, and we assume that the stress 



depends on the instantaneous strain and vorticity. We wish to mention briefly 
that another, non-classical, approach to tides is also possible. In modern books 
on Elasticity (Eringen, 1968) the principle of heredity is re-introduced. It was 
first formulated by Boltzmann (1874) and Volterra (1930). Its basic assumption 
is that the present state of the system, including its stress, depends on all of 
the system's past history. The constitutive equations must then contain a con- 
volution with respect to time. Surely, the principle of heredity falls outside of 
the frame of Analytical Mechanics, but it was recently developed for the sea 
tides by Munk and Cartwright (1966), (Cartwright, 1968) and successfully 
applied to the spectroscopic analysis of tides at Honolulu and on the coast of 
East Britain. Despite this we do not use the principle of heredity in the present 
work, although we recognize that it deserves considerable attention, because 
the Earth might represent a system with memory. 

Basic Notations 

We introduce the following basic notations: 

i, j, k - the basic unit vectors of a rotating system of rectangular 
coordinates with the origin and directions rigidly fixed 
relative to the "initial” (mean) position of the mantle, 

I = ii + j j + kk - the idemfactor, 

r - the position vector of the Earth particle in its "initial" (mean) 
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x, y, z - the rectangular coordinates of P 0 , 

Q - the instantaneous angular speed of rotation of the Earth, 

Q 0 k the constant part of £2 , 

Q 0 IS - the effect of nutation in ft, 

ft = n 0 (k + IN), (1) 

p - the initial density at P Q , 

S - the initial stress at P Q ; we assume that S is a symmetric 
dyadic, 

p - the hydrostatic pressure at P Q , 

r - the initial deviatory stress at P 0 , 

S = - pi + r, (2) 

U(r ) - the initial force function of self gravitation (per unit of mass) 
at IJ ; the force function of the centrifugal force is not 
included, 

R - the position vector of the displaced point P at the moment t, 

w - the absolute acceleration of the displaced point P at the 
moment t, 

u = R - r - the elastic displacement of P 0 at the moment t. 


p' - the density at P in the moment t. 
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V - i _l_ + j _ + k ~ - the del operator relative to r , 

dx oy dz 


e=-(uV + V u ) - the strain dyadic, 

(3) 

oj = y V x u - the local vorticity 

(4) 

ii 

<1 

s 

(5) 

then 


lx co — co x I = -i- (uV - Vu), 

(6) 

(I x co)' = - I X (d 


and we deduce from (3) and (6): 


Vu = 6 - I x w , 

(7) 

uV = e + I x to , 

(8) 


s - the stress increment added to S at the displaced point 
after the deformation caused by the vorticity, 

G - the constant of gravitation, 


- the mean radius of the Earth. 
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Differential Equations of Tidal Oscillations 

Let V be the volume of a portion of matter of the Earth in the initial state 
and A be its boundary. At time t they become V' and A 1 , respectively. The 
corresponding volume and the oriented surface elements we designate by: 

dV, dA, dV‘, dA'. 

Under the influence of the vorticity w a dyadic ab receives the increment 

(cox a)b + a (co x b ) = co x ab - al) x co, 

assuming that « is small. Consequently, under the influence of the vorticity the 
initial stress S receives at time t the increment: 

CO x S - S X CO 

and the total stress becomes: 

T = s + S+ cox S-Sxco. 

The principle of D’Alembert for the elastic tidal oscillations of the Earth in the 
rotating system coordinates takes the form: 


L " 


U(R) + V R (W + -A) 


- w] dV' + J dA' - T = 0, 


(10) 


where 


w =T5.+ 2S2 x f!5 + -^x R + £2 x (£2 X R), 
dt 2 dt dt 


( 11 ) 


is the absolute acceleration of the displaced point, W is the luni solar tidal force 
function, and 0 the force function of the additional attraction, as caused by the 
tidal changes of the internal density and by the formation of tidal "buldges” on 
the surfaces of discontinuity and on the free surface. 



12 


Taking into account: 


BR BR , , BR BR , , BR BR . , 

clA = - — x - — dydz + : — x : — dzdx + - — x - — dxdy, 

By oz Bz Bx Bx dy 


( 12 ) 


and 


dA = idydz +- jdzdx + kdxdy, 


(13) 


we obtain: 


dA' = dA • M, 


(14) 


where the dyadic M has the form: 


M = i 



BR 

Bz 



BR , BR BR 

k aT 37 


(15) 


Substituting 


R = r + u 


into (15) and taking the relations 


I X i - kj - jk , 


I x j = ik - ki 
lx k = ji - ij 

into account, we obtain after some easy transformations: 


M = I - (I x V) x u, 
on, making use of the decomposition 


(I x V) x u = uV - IV • u 


and of (8), we have 


(16) 


M = (1 + 0) I - 6 - I X «. 


( 17 ) 



13 


Making use of the law of conservation of matter 

p' dV* = pdW, 

We can re-write the principle of D'Alembert (10) in the form: 

CVU(r) + u-WU(r) + V(W + 0) - w] dV (18) 

+ J dA- (M*T) = 0, 

neglecting the terms of higher order in u. Both integrals in (18) are now taken 
over the initial volume and surface, respectively. By applying Gauss' theorem 
to (18) we obtain the equation of motion: 

pv = p [V(U + W + + u • VVU] + V* (M* T). (19) 

It is of interest to note that p in (19) is the density corresponding to the initial 
state. From (9) and (17) we have, neglecting terms of higher order, 

M* T = (1 + #)S + s-Sx « - e • S ( 2 °) 

and from this relation there follows: 

V- (M-T) = (1 + 0) V-S - (V-e - V0)- S + V*s 
- (S ’ V) x w - (V * S) x to - € • ’ VS. 

We have from (3) and (4) 

V • e - V0 = + i (V 2 u -W‘u) 

z - - Vx Vx u = — V x w , 

2 


( 22 ) 
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and (21) becomes: 

pw = pV(U -t- W + \p ) + pu • VVU + V • s 

+ (1 + 0) V‘* S + to x (V * S) + (V x co) ■ S 

- (S ■ V) x to - e ■ • VS. 

This last equation is the vectorial form Of the Biot (1965) differential equation 
corresponding to the case of tidal oscillations. 


We can set in (11): 


Q = n 0 (k + N), 


where H 0 is a constant and H 0 N is the effect of nutation. Substituting 

dR = du 
dt dt 

and (24) into (11) and neglecting the insignificant terms, we have 


dt 2 0 dt 0 dt 


+ x (k x r ) + fij|jk x (k x u) 
+ Qgk x (IN x r) + O'- IN x (k x r), 

or, expanding the double products, 


(24) 


(25) 


w - d2u + 20 k x du + 0 dN x r - V<£ 
w--_ + 2U 0 kx _ + iZ 0dt x r V<p 


( 26 ) 


-n%(l - kk ) ’ u + Qq (kN • r +INk'r -2rk-iN) 
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where 


0 = + y 2 )- 


(27) 


There follows from (27): 

W0 = O 2 (I - kk) 

and taking into account 


V [(r x k) • (I\ x r)] = V( r -Nk*r - r 2 k-N) (28) 


= kIN-r + IN k • r - 2rk *N, 
We can re-write (26) in the following form: 


w 


= + 2fi 0 k X ^ - V X 

dt 2 dt 


dN 

dt 


- V[0-flJ(r x k)-(Nx r)] - u • VV0. 


(29) 


Substituting (29) into (19) and setting 

V = U + 0 


We have: 


/7^ + 2Q 0 kx^). :/ oV[V + W + ^-n 2 (rxk)-(Nxr)] + pu • WV 


( 30 ) 




dN 


r x — 


dt 


+ V*s + (l + 0)V*S + tdx (V • S) + (Vx<o)*S - (S’V)x to 


e * - VS. 
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Taking into account the condition of static equilibrium in the initial state 


We obtain from (30): 


V-S + pVV = 0, 


p 





= - p6\J\ +/?(!! • WV ~ w x VV) 


(31) 


+ p V [W + 0 - Qq (r x k) • (N x r )] + V - s 

HN 

+ pQ 0 r x ^ + (Vx co) * S - (S * V) x to - e • • VS. 

We can split the right hand side of the last equation into the main and perturba- 
tive terms by making use of the identity: 

u-WV - W x VV = V(u-W) - e -VV, (32) 

which can be easily proved using (3)-(6). Substituting (32) into (31), we have: 

p (— + 2fi.k x — \ = V-s -pSVV (33) 

\,dt 2 d y 

dIV 

+ pVK + pfi- 0 T X -j-— . — (p£ *W + 6 ■ * VS) 

+ [(V x co) • S - (S • V) x to] , 

where we set for brevity: 

K = W+i/r + u*VV-nJ(r x k)-(I\ x r ). 

The last two terms disappear if the initial stress is reduced to the hydrostatic 


pressure. 
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Together with Biot (1965) we assume a linear relation between the additional 
stress s and the strain e of the form: 

s = B * • £, 

where B is a tensor of the fourth rank with some symmetry properties. Taking 
the second law of thermodynamics into consideration Dahlen (1972) obtains: 


B = C + ~ (IS - SI), 


(35) 


where C is a tensor of the fourth rank, independent of the initial stress. The 
components of C define the set of linear isentropic coefficients. 

Substituting (2) into (35) we have 


B=C +- (It-tI). 


(36) 


If we consider the Earth to be elastic and isotropic, the number of isentropic 
coefficients is reduced to two standard ones, X( r ) and /x(r ). In this case we 
have 


C = XI I + 2yu(i Ii + j I j + k Ik). 
From (34), (36) and (37) we obtain: 


(37) 


s = cr + - (Ir . . € - 9r), 


(38) 


where cr i s the standard stress dyadic, 


cr =■ XI 6 \ 2 jj,e. 


(39) 
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The second term in (38) represents a correction for the deviatory prestress. 
Substituting (38) into (33) and making use of (2), we obtain: 


P 



+ 20 Q k x 



- V-cr - p#VV + pV K 


(40) 


+ pQ 0 r x - e . (pVV - Vp) + D, 

where 

D = A [V (r . . e) - V(6t)} + (V x u) ■ r - (r • V) x co - e . . Vr. (41) 

represents the ’'cross-effects" between the strain and the deviatory pre-stress. 
The last two terms in (40) are perturbative and disappear together with the 
initial deviatory stress. Thus, if we assume that the initial state is in hydro- 
static equilibrium we have simply: 


— + 20, k X — = -v-cr - 0VV + VK + Or x (42) 

dt 2 0 dt p * dt 

This last equation is a slightly modified form of the basic equation of the 
Molodensky tidal theory. The differential equation (40) serves as a foundation 
for the extension of Molodensky theory to the case when the Earth is laterally 
heterogeneous and the initial equilibrium is not a hydrostatic one. 

Similarly, again making use of (2) and (38) and the identities 

pu * VW = V (pu • VV) - (u • VV) Vp - p e ■ VV + pw x VV, 

p0 - V- (pu) - u • Vp, 
u x (VV x Vp) = VV (u • Vp) - Vp (u * VV), 
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we reduce equation (31) to the form; 


/ d 2 u on , du 

W * ^ 


= - V-(pu) VV + V(pu-VV) 


(43) 


dN 


+ pV^ + pQ 0 r x — + V-cr 


V- 


+ u x (VV x Vp) - e • (pVV - Vp) + D, 

where 

V=W + ^- n 0 (rx k)'(Nx r). 


If the initial equilibrium is hydrostatic then the perturbative terms disappear 
and equation (43) becomes: 


d 2 u 

pit 2 


+ 2 n 0 kx^ 


= - V* (pu) VV + V(p U-VV) 


(44) 


4- pW + pQ 0 r x ~ + V ■ cr. 

If, in addition, the geostrophic force and the effects of nutation are neglected 
we have : 

p = - V- (pu) W + V (pu • VV) + pVW. ( 45 ) 

dt 2 

This last equation can serve as a foundation for the theories of oscillations of 
the Earth as developed by Takeuchi (1950), Aiterman and al. (1959) and Saito 
(1971). A simplification of (40) and (42)-(44) is possible if we assume, together 
with Molodensky and Melchior, that the polar component of N does not have any 
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influence on the tidal oscillations of the Earth. This takes place, for example, 
(Melchior and Georis, 1968) if the Earth is an ellipsoid of rotation and the in- 
ternal density depends only on the radius -vector and the latitude. In this case 

k *N = 0 

and, as a consequence, 

(rx k) • (IN x r) = r - kN'r, 

= W ngr-kN-r, 

and to be substituted in the differential equations takes the form: 

, W = V(W i \p) - (kIV • r + INk • r ). 

If, in addition, we consider only terms of a given frequency a, where a is a 

linear combination of the basic tidal frequencies, then in accordance with the 
customary procedure we change the notations and replace 

ubyue 11 ', Nby!Ne iat , 

^-by iak x Ne iat , 
d t 

r x — by iar x (k x IN) = 2iakr ’IN - ia.V(r ’klN'r) 
d t 

where the factors at e iat depend upon the position only. Substituting these 
values into the differential equations, we obtain a linear partial differential 
equation satisfied by the new u: 

p (- a 2 I + 2 iaflgk x I) ■ u = V - a - p6W + pVQ (40 T ) 


+ 2iapQ 0 kr *1N - € • (pVV - Vp) + D, 
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where 

Q = W + 4 1 + u *VV - ^ 0 (^ 0 + ia) r *k N ‘r, 
or 

p(- a 2 I + 2iaO Q k x I)-u = V*cr - V* (pu) VV ^ ) 

+ V (/Ou * VV) + pVQ + 2/Oian o kr - IN 
+ u x (VV x Vp) - e • (pVV - Vp) + D. 

With the corresponding simplifications of (42) and (44) we have: 


(- a 2 I + 2 iaO Q k x 


I)- u = — V*cj - O'vV + VQ 
P 


(42 r ) 


+ 2ia.Q 0 kr *N 

for Molodensky theory, and 

(- a 2 I + 2 iaQ Q k x I) • u - [- V’ (pu) VV + V-cr (44 ) 

+ V(pu-VV)]+ VQ -V 2iafi 0 kr *N 

for the extension of the Alterman-Jarosch-Pekeris theory. Multiplying (40') 
and (42') -(44') by the dyadic 

r = I v i/3k x I - /3 2 kk, 


2Q 0 

p = — . 

a 


we can solve with respect to u and obtain: 

- a 2 (1 — yd 2 ) u — r-P. 

where P is the right hand side of one of the equations (40') or (42 f )-(44 f ). 


( 46 ) 
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The differential equation (46), combined with the equation 

VV = + 4ttGV- ( pu), < 47 ) 

can serve as a basis for the numerical integration of the tidal effects, starting 
from the surface and going, "down." A number of well known continuity condi- 
tions must be satisfied in the transition for one layer to the next. 

In addition, the form (46) helps to recognize the role of the divisor 1 - fi 2 
which becomes small for a - 2H 0 . 


On the Exterior Tidal Potential 

The exterior tidal potential (Biot, 1965) acting on the satellite has the form: 

n' = - G f v,( - ou) Hv + >: G f , (48) 

Jv |r - t ' I J s |r - r f | 

where r is the position vector of the satellite and r' is the position vector of the 
Earth particle. The first integral in (48) is taken over the volume of the whole 
tidally undisturbed Earth. The integrals under the summation sign are taken 
oyer the surfaces of discontinuity, and Sp is the jump in density in transition 
from one layer to the other. We can assume, without loss of generality, that u 
is the solution of the basic partial differential equation and that the periodic 
factor e ,at is omitted. 


We have: 


= 477 


m =+ n 

r' n V" 1 


r - r 


y _ j _ 

£ I 2n + 1 r n + l Z I 

n = 0 m— “n 


Y*(S, a.5)Yj', V), 


(49) 
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where O') are the polar coordinates of the Earth particle relative to the 

system connected with the Earth, (r, a, S) are the polar equatorial coordinates 
of the satellite, and 6 is the sidereal time, and are the normalized com- 
plex spherical harmonics. Substituting (49) into (48) we obtain for the exterior 
tidal potential in the system connected with the Earth: 

I! srr. ' £ £vci>.-* 

n a 0 m = ~ n 


where 


K 


rim 



Y n /. V) V- (pu)dv 


(50) 


+ * j[ ^')^u-dA. 

The numbers K are the analogues of Love elastic parameters. The choice of 

nm 

the initial conditions and of the coefficient of the periodic factor associated with 
u can be arranged in such a manner, that K ^ will depend on the degree n 
stronger than on the frequency and the order n. The exact character of this 
dependence can be established only by the numerical integration of the differ- 
ential equations of the tidal oscillations using an Earth model. Such an integra- 
tion should be planned and performed in the near future. 
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Conclusion 

One of our goals in the present work was to show a connection between the 
exterior tidal potential acting on the satellite and the tidal bodily and surface 
displacements of the solid Earth. We also emphasized the role of static and 
tidal variations in internal density, as well as the role of layers of discontinuity 
in the formation of the exterior tidal potential. The parameters of the Earth’s 
elastic response to tidal forces of a given frequency appear in the present ex- 
pansion of the tidal potential in a natural way. They are represented as a sum of 
the two terms. The first term reflects the influence of the tidal change in the in- 
ternal density and the second term carries the influence of the tidal ’’buldges." 
The dependence of the elastic parameters (Love numbers) on the local coordi- 
nates on the Earth surface may, in fact, reflect the asymmetric and inhomo- 
geneous structure of the interior of the Earth.*) The differential equations of 
tidal oscillations we develop in the present article are valid for the laterally 
heterogeneous Earth and for the presence of initial deviatory stress. 

Like ail equations of other theories of tidal oscillations they can be inte- 
grated only numerically. 


*D. E. Smith. Private communication. 1972. 
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